Square vortex solitons with a large angular momentum 



O 

o 



Humberto Michinel, Jose R. Salgueiro and Maria J. Paz-Alonso 

Area de Optica, Facultade de Ciencias de Ourense, 
Umversidade de Vigo, As Lagoas s/n, Ourense, ES-32005 Spain. 

We show the existence of square shaped optical vortices with a large value of the angular mo- 
mentum hosted in finite size laser beams which propagate in nonlinear media with a cubic-quintic 
nonlinearity. The light profiles take the form of rings with sharp boundaries and variable sizes 
depending on the power carried. Our stability analysis shows that these light distributions remain 
stable when propagate, probably for unlimited values of the angular momentum, provided the host- 
ing beam is wide enough. This happens if the peak amplitude approaches a critical value which 
only depends on the nonlinear refractive index of the material. A variational approach allows us 
to calculate the main parameters involved. Our results add extra support to the concept of surface 
tension of light beams that can be considered as a trace of the existence of a liquid of light. 

PACS numbers: 46.65.Jx, 42.65.Tg 
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I. INTRODUCTION 

In wave mechanics, a vortex is a screw phase disloca- 
tion, or defect Q where the amplitude of the field van- 
ishes. The phase around the singularity has an integer 
number of windings, I, which plays the role of an angular 
momentum. For fields with non- vanishing boundary con- 
ditions, this number is a conserved quantity and governs 
the interactions between vortices as if they were endowed 
with electrostatic charges 2J. Thus, I is usually called the 
"topological charge" of the defect. 

Vortices are present in very different branches of 
physics like fluid mechanics, superconductivity, Bose- 
Einstein condensation, astrophysics or laser dynamics, 
among othersQ. In Optics, a vortex with charge I takes 
the form of a black spot surrounded by a light distri- 
bution. Around the dark hole, the phase varies from 
zero to 2irl. These defects appear spontaneously in light 
propagation through turbulent media and can also be 
produced by appropriately shining a computer generated 
hologram|j|. The trace of vortices in a light field is a 
characteristic "fork-pattern" interferogram produced by 
superposition with a tilted planar wave. 

The first experimental works on optical wavefront dis- 
locations were carried out in the 80's, in the context of 
adaptive systems, where phase singularities were a severe 
problem for image reconstruction techniques^. Later on 
they have been studied, among other fields, in optical 
tweezing6], particle trapping 7], laser cavities 8], opti- 
cal interconnectorsj^ or even to perform N-bit quantum 
computers [Toj. 

Concerning light vortices in the nonlinear regime [TT| . 
the first theoretical work analysed their stability in 
Gaussian-like distributions propagating in optical Kerr 
materials|l2j. It was found for a cubic self-focusing re- 
fractive index, that a beam of finite size will always fila- 
ment under the action of a phase dislocation. This also 
applies to saturable self- focusing nonlinearities [l3j . On 
the other hand, vortex states were predicted and found 
experimentally for self-defocussing materials both in the 
Kerr case for continuous background^) and in the sat- 



urable case with finite size beams [15J. 

It was shown in 16] that stable vortex states with 
I = 1 can be obtained as stationary states of the prop- 
agation of a laser beam through cubic-quintic optical 
materials ^3 0, [^. This kind of nonlinearity is char- 
acterized by the > and < components of 
the nonlinear optical susceptibility and changes from self- 
focusing to self-defocusing at a given intensity [2(j ■ It has 
been recently shown that a gas-liquid phase transition 
takes place in light beams propagating in this type of 
materials pi). 

In this work, we will show that stable vortex states 
with a huge value of the angular momentum exist and 
their peak amplitude and propagation constant tend 
asymptotically, as the beam flux is increased, to val- 
ues that do not depend on In this way, our results 
are in contradiction with previous workp^]. where it was 
claimed that stable vortex states in finite size beams ex- 
ist only for the values I = 1,2. For I = 3 it was found 
a persistent weak instability which was also supposed to 
exist for higher values of the angular momentum. 

In next section we will analyze the cubic-quintic non- 
linear model, finding numerically the stationary states for 
a wide range of the angular momentum / (up to 50) and 
describing their particular properties. Then, we will cal- 
culate analytically, by means of the variational method, 
the critical values of the propagation constant and peak 
amplitude that characterize the domain of existence of 
vortices. Finally, we will perform an azimuthal stabil- 
ity analysis to determine the domain zone where stable 
states can be found. 



II. THE MODEL 

Let us start by writing the equation for laser beam 
propagation along z in an optical cubic-quintic material. 
For paraxial propagation, the equation for the beam en- 
velope ^ is a nonlinear Schrodinger equation (NLSE) of 
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FIG. 1: Numerically calculated radial amplitude profiles of 
the stationary states of Eq. JHJ for I = 3, 4, 10 and 50 with 
/3//3 cr =0.1, 0.5, 0.9 and 0.99. In all the cases n 2 = n 4 = 1 . 



the form: 



.9* 



J— + Vi* + (n 2 |<F| 2 - n 4 |<F| 4 ) * - 0, 



(1) 



where 



2 d 2 /d<t> 2 



d 2 /dr 2 is the 



f^d/dr 

transverse Laplacian operator in cylindrical coordinates 
(r, 0, z). The real positive constants n 2 and ro 4 are given 
respectively by the x'" 3 '' > and < components of 
the nonlinear optical susceptibility and characterize the 
dependence of the refractive index on the intensity of the 
beam. If n 4 < 0, a Gaussian beam of high enough power 
will undergo collapse after seli-focusing|23j. The effect of 
a negative fifth order susceptibility (— n 4 term) combined 
with diffraction will stop the collapsing tendency for high 
powers, yielding to a stable two-dimensional condensed 
state of light with surface tension properties similar to 
those of usual liquids HJ E H3- 

We are interested in stationary states with radial sym- 
metry and angular momentum I of the form: 



(2) 



where [3 is the nonlinear phase shift or propagation con- 
stant and tp(r) is the radial envelope of the field. After 
substitution of @ in (JIJ, the following z-independent 
equation is obtained for ip(r): 



I 2 

-PiP + Vty ~^ip + n 2 ip 3 - n 4 V 5 = 0, (3) 

where V 2 = d 2 /dr 2 + (l/r)d/dr, is the radial part of the 
Laplace operator. 

For a given integer value of I, a continuum of eigen- 
states with ip — > as r — ► oo can be obtained by solving 
numerically Eq.©. Close to the origin the shapes fol- 
low the linear regime with ip oc r l . To this aim, we have 



used a standard relaxation technique. The profiles of the 
eigenstates for several values of I and (3 are plotted in 
Fig. \I\ior the case of n 2 — n 4 = 1. We particularly show 
states with I — 3 and I — 4 since these were previously 
found unstable in previous work [22j, as well as two ex- 
amples of large angular momentum states (I = 10 and 
I = 50). In all cases, the stationary states can only be 
found for values of /3 between zero and a fixed critical 
value /3 cr 0|, which does not depend on I. 

It can be appreciated in the graphs that values of (3 
below 0.5/3 cr yield to light distributions with smooth and 
wide Gaussian-like shapes. As (3 is incremented the beam 
flux grows and the spatial profiles narrow, yielding to a 
minimum thickness of the ring of the stationary state 
for values of (3 around 0.8/3 cr , keeping approximately the 
Gaussian shape. For larger values of the propagation 
constant, the beam flux grows rapidly with (3 and the 
peak amplitude of the distribution saturates due to the 
effect of n 4 , reaching asymptotically the value A cr which 
is slightly below the maximum amplitude. Thus, high 
power beams show spatial light distributions with flatted 
tops in their profiles, similar to those of hyper-Gaussian 
functions IMISH, 
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FIG. 2: Maximum amplitude of the stationary states versus 
propagation constant for angular momenta 1=0, 1, 2, 3, 10 
and 50. All the curves join at the point A = A cr = 0.866, 
j3 = f5 cr = 0.1875. Inset: detail of the zone close to /3 cr where 
the calculations become delicate. 

We must stress the intriguing fact that both j3 cr and 
A cr do not depend on the value of the topological charge. 
This is shown in Fig. [2 where the maximum amplitude 
has been plotted as a function of (3. In the inset, the zone 
[3 w [3 C r can be seen in detail. As can be appreciated, 
whatever the value of I is, all the curves tend to join at 
the same point. This means that the critical value of the 
propagation constant and peak amplitude only depend 
on the nonlinearity and not on the angular momentum. 
We will revise this result in our analytical study of the 
next section. 

It also worths to mention that the central hole increases 
its size with the topological charge for a fixed value of (3, 
as can be seen comparing the profiles in Fig. I = 3, 4 
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FIG. 3: (Color on line) Azimuthal eigenstates of Eq.J^Jl for 
I = 1 to 9 with = 0.95/3 cr . 



with I = 10, 50. This is also clearly shown in FigEl where 
we plot several eigenstates with values of the angular mo- 
mentum ranging from / = 1 up to I = 9, with propagation 
constant (3 — 0.95/3 cr . Besides, if (3 grows, the radius of 
the hole increases. As the value of (3 approaches j3 cr the 
thickness of the external ring grows faster than the inter- 
nal hole, and the final result takes the asymptotic form of 
a dark spot surrounded by a larger ring of light of almost 
constant shape which ends abruptly at a given radius. 
This behavior can be assessed having a look to FigQ] 
where the dimensions of the internal hole and the ring 
thickness are plotted versus the propagation constant for 
the particular case of / = 10. A logarithmic scale was 
chosen to highlight that the growth in the ring thick- 
ness clearly dominates over the hole radius from certain 
value of the propagation constant. In the inset it is also 
shown, as an example, one of the stationary states with 
(3 very close to (3 cr , showing the huge ring whose width 
clearly exceeds the hole radius and presents a practically 
rectangular shape. 



III. VARIATIONAL ANALYSIS 

To explain the properties of the above light distribu- 
tions, we have performed a variational analysis pi l25|. 
It is easy to demonstrate that the stationary system de- 
scribed by the model © can be obtained from the fol- 
lowing Lagrangian function: 
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FIG. 4: Dependence of the internal hole radius (continuous 
line) and ring thickness (dashed line) with the propagation 
constant j3 in the vicinity of (3 cr for eigenstates with I = 10. 
Inset: example of a stationary state with / = 10 and /3 = 
0.9998,3 cr .. 



Now assuming a rectangular shape of the stationary 
states for values of (3 close to (3 cr (see Fig. ^ , we choose 
a trial function ip v (r) centered at ro with amplitude A 
and width 2w, given by: 



i>v(r) 



A, |r — ro| < w 
0, |r — r | > w 



(5) 



Thus, ro, A and w are the variational parameters. Using 
this trial function and minimizing the average over the 
Lagrangian < C > respect to each parameter |26j. we 
obtain the following conditions: 
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13 = 0, 
-f3 = 0. 



(6) 
(7) 



where Eq. JBJ is obtained from the minimization respect 
to ro and w (the same condition is obtained for both 
parameters), and Eq. Q follows from minimization re- 
spect to parameter A. In the limit when (3 — > /3 cr , the 
first term of both equations vanishes. This follows by 
taking into account that r > w should always be sat- 
isfied (otherwise there would be no hole), then we have 



-w > ri 



oo and (r. 



0, consequently the 



first term of Eq. © is zero. For Eq. J7J, the argument of 
the logarithm tends to infinity as it is easily deduced from 
the fact that the ring width grows faster than the hole ra- 
dius and ro > w ((ro + w)/(ro — w) > 2w/(ro — w) — > oo). 
However, this term diverges logarithmically, meanwhile 
the denominator goes to infinity quadratically (product 
row), and consequently the whole term tends to zero. 
Finally, we can solve them for A cr and p cr to obtain: 
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For the particular case considered in the numerical cal- 
culations displayed in Figs. 11121 i.e. taken n 2 = ri4 = 1, 
the values obtained for the critical parameters are A cr = 
0.866025 and (3 cr = 0.1875. The comparison of these an- 
alytical results with the numerical calculations shows an 
excellent agreement, since both values are exactly those 
guessed numerically (see Fig- EJl - This so good result is 
due to the choice of the trial function, which fits almost 
exactly with the numerical solution for values close to the 
critical point. 



IV. STABILITY ANALYSIS 

In order to test the stability of the stationary states we 
calculated the growth rates of small azimuthal perturba- 
tions to find out the value of f3 at which they vanish. 
Additionally, in order to assess the accuracy of the pre- 
vious analysis, we propagated some unstable eigenstates 
with a split-step Fourier method and found their split- 
ting distances. The inverse of these values should coin- 
cide, except for a constant scale factor, with the dominant 
perturbation eigenvalues calculated in the azimuthal in- 
stability analysis. Finally, we have also simulated other 
kind of perturbations like total reflection at the bound- 
ary between a cubic-quintic material and air. As we will 
see below, the eigenstates show robust behavior against 
these collisions and preserve their angular momentum al- 
though strong oscillations are observed. 

To carry out the perturbation analysis, we add to the 
original eigenstate a small p-order azimuthal perturba- 
tion function |13ll27j : 



§(r, 0, z) = [ip(r) + f(r, z)e^ + h(r, z)e^\e li - l ^ Siz \ 

(10) 

where /(r, z) and h{r, z) are the small complex compo- 
nents of the eigenstate of the p-order azimuthal perturba- 
tion. Our interest is to seek those functions which grow 
exponentially with z, so we assume that they have the 
form: 



f(r,z) = [f 1 (r)+if 2 (r)}e 5 **, 
h(r,z) = [h^r) + ih 2 (r)]e- s > , 



(11) 
(12) 



being the parameter 8 P the perturbation eigenvalue. In 
this way, the real part of 6 P constitutes the growth rate 
of this perturbation. If we replace the perturbed eigen- 
state (Eq. COl) into Eq. Q and keep only the first order 
terms in /(r, z) and h(r, z) (linearisation) we obtain the 
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FIG. 5: Growth rate of the azimuthal perturbation versus (3. 
In (a) and (b), the angular momentum is fixed (I = 3 for (a), 
1 = 4 for (b)) and the values of the perturbation order p are 
indicated by labels near the curves. In (c) the growth rate 
for the perturbations with p — 2 for I = 1, 2, 3, 4, 10 and 50 
are plotted. In (d), a comparison between the inverse of the 
splitting distance (dots) and p = 2 perturbation eigenvalue 
for unstable states with I = 3 is shown. 



following set of coupled differential equations for those 
components f{r,z) and h(r,z): 



# + V?/-^/ + W)/ + ]J(#* = (13) 

OZ T A 



i— +V 2 r h 

OZ 



Q—PL h + Q(il>)h + RW)r=0, (14) 



where Q(V>) = -/? + (2n 2 - 3n 4 \^\ 2 )\ip\ 2 and R(ip) = 
(ri2 — 2tt,4|'0| 2 )|V'| 2 . The solution of this equation system 
is obtained using a Crank-Nicholson scheme to propagate 
an initial arbitrary guess until the shape of each compo- 
nent does not change perceptibly 27]. According to the 
component dependence on z (Eq. [fl"T|l h the value of the 
growth rate can be calculated at each propagation step 
by: 



1 , \f(r,z + Az) 



In 



(15) 



where Az is the propagation step and the function /(r, z) 
is evaluated in a fixed point r, usually that which corre- 
sponds to the maximum. Besides, the functions can be 
rescaled at each step by this maximum value to avoid an 
overflow. The propagation is carried out until the value of 
the perturbation growth rate does not change any more, 
what indicates that convergence was reached. This al- 
lows us to obtain the growth rates Re(r5 p ) for different 
order perturbations versus the propagation constant, as 
depicted in Fig. |SJ 

The growth rates for vortices with angular momen- 
tum I = 3 and I = 4 are shown in Figs. Ela)-(b). As 
can be seen, all of them fall to zero for a value of (3 
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FIG. 6: (Color on line) Numerical simulation of the total 
reflection at a planar boundary between a cubic-quintic and a 
linear media (air). For this stationary state the angular mo- 
mentum is I — 4 and the propagation constant /3 = 0.95/3 cr . 
The top image is an isosurface of the propagating beam. The 
boundary between the nonlinear material and air is the plane 
y — 0. The images (a)-(c) correspond to different transversal 
cutoffs (intensity profiles) of the beam at different values of 
the propagation distance z. 



below (3 cr , what implies the existence of a stability win- 
dow, in contradiction with previous calculations where 
all the states with I > 2 where found unstable j^. Our 
results show that the maximum growth rate corresponds 
to perturbation eigenvalues with p m 21, what allows to 
estimate the number of filaments N resulting from the 
breakup of the unstable vortices (N « 21). Besides this, 
the perturbation p = 2 has been proved to be the most 
persistent, despite the value of the angular momentum. 
Hence, in Fig. [SJc) we plot the curves associated to this 
perturbation for different values of the angular momen- 
tum, including the cases corresponding to I = 10 and 
/ = 50. As can be appreciated in these plots, there ex- 
ists a window between the vanishing point and the limit 
value for (3 (f3 = (3 C r), what proves the existence of a 
stability zone close to the critical point containing an 
inhnite number of stable eigenstates. Note that this win- 
dow narrows for high values of I but remains finite. As 
I increases, the point at which the perturbation vanishes 
approaches asymptotically the critical point. However, 
we believe that the critical point itself is never reached, 
even for arbitrarily high values of the topological charge. 
When (3 is close to j3 cr the azimuthal analysis turns 



itself very delicate and it has to be carried out in a very 
careful way. In fact, convergence takes a much longer 
distance and an erroneous final result is obtained if the 
number of samples and the propagation step are not cho- 
sen appropriately. In this sense, combining the analysis 
with direct calculations of the splitting distance of the 
unstable eigenstates is definitively useful. In Fig. OJd) 
it is zoomed the region oflSJc) where the perturbation 
for I — 3 drops to zero. The points obtained propagat- 
ing the eigenstates and taking the inverse of the distance 
where they split are also plotted. These values were sub- 
sequently scaled by the same constant value to compare 
with the perturbation eigenvalue curve. As can be appre- 
ciated, the values obtained from this propagation experi- 
ments fall to zero with the same slope as the perturbation 
eigenvalues do. When the stability analysis is not per- 
formed with enough accuracy a more steady behavior of 
the curve appears, what implies that the eigenvalue falls 
to zero at a higher value of (3. This allows us to assess 
the validity of the perturbation analysis. 

As a final test of the stability of the eigenstates, we 
have simulated the total reflection at a planar boundary 
between a cubic-quintic material and air for beams with 
different angular momenta. For the simulation we have 
used a split-step Fourier method with a 520 x 520 grid. 
The idea is similar to the test of surface tension properties 
of " liquid light beams" from ref . . As can be seen in 
Fig. El a beam with / = 4 does not split after the total 
reflection, although a strong oscillation is observed. This 
is another proof of the stability of these nonlinear waves. 
We must notice that depending on the incidence angle, 
a strong deformation of the beam can be induced, which 
can yield to a split or a decay of the inner vortex into 
several defects with lower topological charges • 



V. CONCLUSIONS 

The main conclusions that can be derived from the 
present work are the following: first, stable azimuthal 
finite-size beams with arbitrary very large angular mo- 
mentum can exist in optical materials with self-focusing 
(cubic) and defocussing (quintic) nonlinearity. Second, 
the shapes of these beams tend asymptotically to square- 
like ring profiles with bigger dark holes for higher values 
of the angular momentum. And finally, the critical val- 
ues of the propagation constant and amplitude do not 
depend on the angular momentum of the beam. 
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